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Motivated by the search for best lattice sphere packings in Euclidean spaces of large dimensions we study 
randomly generated perfect lattices in moderately large dimensions (up to d = 19 included). Perfect lattices are 
relevant in the solution of the problem of lattice sphere packing, because the best lattice packing is a perfect lat- 
tice and because they can be generated easily by an algorithm. Their number however grows super-exponentially 
with the dimension so to get an idea of their properties we propose to study a randomized version of the algo- 
rithm and to define a random ensemble with an effective temperature in a way reminiscent of a Montecarlo 
simulation. We therefore study the distribution of packing fractions and kissing numbers of these ensembles 
and show how as the temperature is decreased the best know packers are easily recovered. We find that, even 
at infinite temperature, the typical perfect lattices are considerably denser than known families (like Ad and 
Dd) and we propose two hypotheses between which we cannot distinguish in this paper: one in which they 
improve Minkowsky's bound ~ 2~('' *'**'' ''^)'', and a competitor, in which their packing fraction decreases 
super-exponentially, namely <f) ~ d""'^ but with a very small coefficient a — 0.06 ± 0.04. We also find prop- 
erties of the random walk which are suggestive of a glassy system already for moderately small dimensions. 
We also analyze local structure of network of perfect lattices conjecturing that this is a scale-free network in all 
dimensions with constant scaling exponent 2.6 ± 0.1. 



I. INTRODUCTION 

Sphere packing is a classic problem with many connections 
to pure and applied mathematics (number theory and geom- 
etry'), communication theory ' and physics . The statement 
of the problem is very simple: given an Euclidian space of 
dimension d what is the densest spatial aiTangement of im- 
penetrable spheres? In a more formal way one seeks to find a 
maximum over all packings 

</'best(rf) = max(/)(7'). 

■Pes 

Here P is a packing of spheres (an allowed configuration of 
the impenetrable spheres), S is the set of all packings, and 
(f){V) is the fraction of space covered by the packing V. 

As is often the case with problems related to number the- 
ory, the simplest questions do not have simple answers. De- 
spite over 200 years of research the problem has only been 
solved for d = 2 and d = 3 (the famous Kepler's conjec- 
ture). The latter case has only been proven about fifteen years 
ago and required substantial amount of computer work. Al- 
though good and very good candidates for the best packings 
have been identified in higher dimensions (namely < 30) our 
knowledge deteriorates quickly as dimensions become really 
high, say of order 10'^ where the problem becomes of interest 
to communications theory. 

One the greatest challenges in the sphere packing problem 
is that no universal behavior is identifiable. Every dimension 
seems to be peculiar, with some dimensions being very spe- 
cial, like 8, 12, 24. In the generic case there is no restriction 
on packings: they can be of any nature, ordered (crystalline 
breaking of translational symmetry) or even disordered. For 
relatively low dimensions, d < 9 the best (known) packings 
are all lattice packings, that is packings where spheres are 
placed at the vertices of a certain Bravais lattice (one parti- 
cle per unit cell of the lattice). In d = 10 for the first time. 



the best known packing is generated by a non-Bravais lattice . 
Some recent works • conjecture that in high enough dimen- 
sions completely disordered packings might win over regular 
ones. 

To understand the degree of difficulty of the problem it is 
sufficient to mention that even finding good upper bounds on 
best packing fractions uniformly valid for all dimensions re- 
sisted to all attacks so far The one-hundred year old lower 
bound by Minkowsky only received linear improvements un- 
til today and an exponential improvement • only exists sub- 
ject to an interesting but very strong conjecture . Even worse, 
Minkowsky's bound is non-constructive, and no methods are 
known which would allow to construct a lattice which sat- 
isfies at least that bound in very high dimensions. Arguably 
the most important recent contribution in this respect has been 
given by the works ■ in which the problem is reduced, for 
any given dimension, to an infinite linear programming prob- 
lem. The technique is powerful -in 8 and 24 dimensions 
the bounds are saturated by the best known packing, proving 
hence their global optimality- but has not yield an understand- 
ing of the problem for generic d. 

Given the complexity of the generic case it might prove 
useful to consider a simpler version of the problem. One of 
them is the so called lattice packing problem, which restricts 
allowed packings to Bravais lattice packing only. Although 
the set of possible packings is severily reduced, exact results 
are only established up to d < 8, with = 9 case hopefully, 
closed in 2012. 

In theory the lattice sphere packing problem is simpler, be- 
cause it admits an explicit algorithmic solution where one 
has to check a finite number of special lattices to find the best 
one. The best packing, in fact, is both a perfect and eutactic 
lattice (we give the characterization of these lattices later) and 
both the number of perfect ' and that of eutactic lattices is 
finite (hence the intersection is). This algorithm has been 
applied to dimensions d < 8 to systematically find all such 
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lattices"" In this paper we will ran a randomized version 
of the algorithm in dimensions 8 to 19 to generate large (up to 
several millions) set of perfect lattices in each dimensions and 
then study the statistical properties thereof. We will introduce 
a fictitious temperature to explore non-typical regions of the 
space of perfect lattices and get the best known packings. 



II. LATTICES, PERFECT LATTICES AND EUTACTIC 
LATTICES 



A. Notation 

In this paper we will consider only lattices or in Physics 
terminology Bravais lattices, namely lattices which have only 
one particle per unit cell. A generalization of our results to 
arbitrary but finite number of particles per unit cell will be 
discussed at the end of the paper. In our definitions and logic 
of discussion we will follow closely Schurmann although 
we will not pretend to achieve the same level of rigor. 

We will define a lattice A, one particle per unit cell, in M'* 
by means the square matrix of the components of the d, d- 
dimensional Unearly independent (basis) real vectors e* 



(e\ e\ ... e\\ 



(1) 



\ei ei ... e^J 
The points in the lattice are elements of the set 

A = {x : X = Az, z e Z'^/{0}}. 



(2) 



The associated symmetric, positive definite d-hy-d quadratic 
form Q is defined by matrix multiplication as 



Q = A^A. 



(3) 



We will refer without difference to the quadratic form Q or to 
the basis matrix A when we talk about a lattice. The distance 
of a point A-z in the lattice is (here T stands for transpose, both 
of a vector and of a matrix) 



I 



Vz'T ATA z = Vz^Qz, 



(4) 



where z^Qz = Y.i.j=i ^iQijZj. 

The notion of shortest vector of a lattice is fundamental in 
the theory of lattices and allows one to connect to the theory 
of sphere packing. Namely define the arithmetic minimum of 
a lattice Q as square of the minimum length of a vector in the 
lattice 



A(g) 



mm 

zeZ'*/{o} 



z^Qz, 



(5) 



and the set 



Min(Q) = {z e Z'^ : z^Qz = \{Q)) 



Let us point out that the set Min(Q) should contain at least 
two vectors (as x and — x have the same length) but for the 



"interesting" lattices the cardinality of the set (known as the 
kissing number) is usually much larger, sometimes even ex- 
ponential in d. The maximum cardinality of Min((5) over the 
set of d-dimensional lattices is an open problem in most d and 
has been dubbed the kissing number problem . 

The connection with the sphere packing problem is eas- 
ily made. The largest non-overlapping spheres we can fit in 
a lattice must have as radius half the length of the shortest 
vectors of Q. Considering that the volume of a unit cell is 
det A — -\/det Q we have the maximum fraction of space 
covered by a sphere packing Q is the ratio of the volume of 
this sphere divided by the volume of the unit cell: 



det(Q)i/2 ■ 



where Bd is the volume of a d-dimensional unit sphere 



B, 



dT{d/2)' 



(7) 



(8) 



A strictly related quantity is the Hermite constant of Q (in 
terms of which the packing fraction can be expressed) 



H{Q) 



HQ) 



(9) 



In the following we will also use another indicator that we 
will call "energy" as a target function to minimize with the 
introduction of a temperature: 



e(Q) = -ilog(</)(g)). 



(10) 



Minkowksy's bound ensures that this quantity is bounded on 
the best lattices even in the limit d ^ oo. 

The lattice sphere packing problem (henceforth LSP prob- 
lem) in d dimensions is the problem of finding the maximum 
of 4>{Q) (or H{Q)) among all the d-dimensional lattices. The 
problem is solved for d = 1, ...,8 " and d = 24 ' only. 



B. Perfect lattices 

We will now concentrate on a subset of lattices which turns 
out to be fundamental in the solution of the lattice sphere 
packing problem: the perfect lattices. 

A lattice is named perfect iff the projectors built with its 
shortest vectors span the space of symmetric d-hy-d matrices. 
So for a perfect lattice Q let Z be the cardinality of Mm{Q) 
and let Vq G Mm{Q), a — 1, Z (Z is also called the kiss- 
ing number of a lattice). Let M be any symmetric d-hy-d 
matrix there exist a set of reals pa such that: 



(6) For example take the square lattice in d = 2: 
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1 



(11) 
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the shortest vectors are 

Min(Q,,) = {(l,0),(0,l)} 
and the projectors are 



Pi = 



1 




P2 = 




1 



(13) 



(14) 



which do not span the space of symmetric matrices. Therefore 
the square lattice is not a perfect lattice. 
Instead, consider the hexagonal lattice 



Qhex — 



2 1 
1 2 



(15) 



It has three shortest vectors (of length \/2) 



MiniQhe.) = {(1, 0), (0, 1), (1,-1)} (16) 
and the corresponding projectors are 



Pi = 



1 




P2 = 




1 



P:^ = 



-1 1 



(17) 

and the reader can verify that they form a basis for symmetric 
2-by-2 matrices (one can easily form linear combinations of 
the P's to obtain the identity and two of the three Pauli ma- 
trices). Note that the number of shortest vectors of a perfect 
lattice is bounded from below by {d + l)d since this is twice 
the smallest possible number of projectors that can span the 
space of symmetric matrices (the dimension of space of sym- 
metric matrices). So in the previous example we could have 
said beforehand that the square lattice is not perfect but we 
should have checked anyway that the hexagonal lattice was 
indeed perfect."^ 

Voronoi proved' ^"'^ that perfect forms are vertices of the 
Ryshkov polyhedron defined as a set of forms Q whose 
shortest vector is larger than a given value: 



Vx^{Q: A(Q) > A} 



(18) 



where the actual value of A (as far as A > 0) is immaterial as 
the axis can be rescaled freely. Therefore we can reduce the 
sphere packing problem on V\, hence constraining to forms 
with X{Q) = A without any loss by finding 



H = 



A 



inf 



(19) 



The number of vertices of the Ryshkov polyhedron, and hence 
of perfect forms is (up to isometries that we define below) 
finite (a small subset of all the lattices in any given dimension 
d). 

The main result which gives importance to perfect lattices 
in the context of the LSP problem is the classic Voronoi 's the- 
orem which can be stated as follows: 



Theorem: 

tice. 



the best lattice sphere packing is a perfect lat- 



The proof (which we do not give here, see'"') follows if one 
shows that det^^'^((5) does not have stationary points inside 



FIG. 1. Left. Square lattice. Right. Hexagonal lattice (also known 
as triangular lattice). 



the Ryshkov polyhedron. This in fact implies that the mini- 
mum of det((3) and the maximum of cj) (or H) occur on the 
vertices of the polyhedron, hence on perfect lattices. 

Therefore the problem of LSP is reduced to finding all the 
perfect lattices and comparing their packing fractions: it be- 
comes a problem for a computer to solve. Unfortunately (or 
maybe, fortunately) things are not so easy as they might seem. 
Indeed the number of perfect lattices grows very fast with the 
dimension (probably faster than exponential, as we will argue 
later) and the task of finding them all has been completed up 
to d = 8 (where they are 10916). For d = 9 has found 5 • 10^ 
forms but the conjectured total number should be around 
2 • 10^. 



C. Isometry of lattices 

A lattice admits many equivalent representations in terms 
of quadratic forms Q: one can rotate the lattice or replace its 
basis vectors with their independent linear combinations. This 
equivalence is captured by notion of isometry: 

Definition: Lattices Q and Q' are isometric if there exists 
a matrix U € GLj^CZ) and c G M such that 

Q' = cU*QU. 



Another name in use is arithmetical equivalence. For example 
the hexagonal lattice Qhex given by Eq. (15) has an equivalent 
representation 

o' -( ^ 

which is isometric to Qhex with isometry matrix 

^ \-i 

A practical way of checking if a given pair of forms are 
isometric was developed in where one uses backtrack search 
to construct an isometry matrix (if this exists). However most 
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of the times it is sufficient to check if some criteria (like the 
number of shortest vectors) are satisfied before running the 
generic code which can be quite slow in high dimensions. 

D. Eutaxy 

The last concept that we need for our investigation is that of 
eutactic lattice. This is not strictly necessary for understand- 
ing our results in this paper but it gives a suggestive connec- 
tion with the theory of spin glasses which we plan to investi- 
gate as a continuation of this work. Eutactic lattices cannot be 
improved (as we will prove below) by an infinitesimal trans- 
formation of the matrix base and therefore are local maxima 
of the packing fraction. Their number also grows with the di- 
mension d and one is then led to think that in high enough 
dimensions this phenomenon is reminiscent of the landscape 
of a mean-field spin glass free energy . 

Given a perfect form Q we can always write (since it is a 
symmetric, nonsingular matrix) its inverse Q^^ in terms of 
the projectors built on its shortest vectors 

Q-^ = J2 ^20) 

xSMin(Q) 

(here xx^ is the matrix with elements XiXj). 

Definition: A eutactic form is one for which one can 
choose all the above > 0. An equivalent definition is that 
is in the interior of the Voronoi domain of the perfect 
form Q, defined as 

V(Q) = cone{xx^ : x e Min(Q)}, (21) 

the cone in the space of forms generated by the projectors built 
with the shortest vectors of Q. 

The Hermite constant (or packing fraction) of an eutactic 
form can only be decreased by any infinitesimal change of the 
form. In fact, by using the identity 

Tr ((V det Q)A) = dct(0)Tr (Q^^A) (22) 

we obtain, to first order in SQ = Q' — Q where Q' e 'Px^q) 
(so the length of the minimal vectors is unchanged): 

H{Q + 6Q) = H{Q) - -^/^Tr {Q-^SQ) < H{Q) 
dct ' (Q) 

(23) 

where the inequality follows from: 

Tx{Q-\Q' -Q) ^ ax(x^O'x-x^Qx) >0, 

xeMin(Q) 

(24) 

as Q' e Vx^Q) and ctx > 0. 

It follows then that a perfect and eutactic lattice is a local 
maximum of H from which 

Tlieorem: perfect and eutactic (PE) lattices are local max- 
ima of the Hermite constant and hence of the packing fraction 
and therefore 



Corollary: the best packing lattice is both perfect and eu- 
tactic. 

The concept of eutaxy is extended to arbitrary lattices with 
introduction of weakly-eutactic, semi-eutactic and strongly- 
eutactic lattices. Weakly-eutactic lattices satisfy Eq. (20) with 
real coefficients a^, semi-eutactic lattices have > (i.e. 
some of the coefficients in Eq. (20) are zero) and finally 
strongly-eutactic lattices are eutactic lattices with all Ux equal. 
Recall that by definition a perfect lattice is (at least) weakly- 
eutactic since xx-^ span the space. The interest in strongly- 
eutactic lattices comes from the fact they are also the best 
packers locally among lattices with arbitrary number of parti- 
cles per unit cell . 

The problem of determining eutaxy class of a form admits 
an efficient solution: given a form, its eutaxy class - non- 
eutactic, weakly-eutactic, semi-eutactic or strongly-eutactic - 
can be decided by solving a sequence of linear programs-- 
and therefore is of polynomial complexity with respect to the 
number of shortest vectors (which, however can grow as fast 
as an exponential of d). 

Summarizing, the take home messages of this section are 
that the maximum of the packing fraction over lattices in any 
given dimension is attained by one of the PE lattices, of which 
there is a finite number (in any given d) and that each of the 
PE lattices is a local maximum. This characterization is ex- 
tremely powerful but still does not prevent us from having to 
find all perfect lattices and checking which ones are eutac- 
tic and which are not. There is a simple and efficient way to 
generate perfect lattices but there is not (as far as we know) a 
similarly efficient way to generate eutactic ' or PE lattices. 
One should first generate perfect lattices and then check them 
for eutaxy. The simple and efficient way to generate perfect 
lattices is given by Voronoi algorithm, which we review in the 
following section. 



ni. VORONOI'S ALGORITHM AND ITS 
RANDOMIZATION 

We have now reduced the problem of finding the best lat- 
tice packing to that of finding the best lattice packing among 
perfect and eutactic lattices. We need a way to generate all 
the perfect lattices, select the eutactic ones and look at the 
most dense among them. The first task is accomplished by the 
Voronoi algorithm • ■ that we now describe. 

Start with a perfect form Q. 

1. Find all the shortest vectors x g Min(Q), and the in- 
equalities describing the cone 7^(Q) 

ViQ) = {Q'\ Vx e Min(Q) : x^Q'x > 0} (25) 

2. Find all the extreme rays of the polyhedral cone V{Q). 
Call them Ri, ...,Rk- 

3. Create the forms Qi — Q + ctiRi, choosing rational 
numbers a; such that the new form Qi is again perfect. 
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4. Check for isometrics and repeat from Start with each of 
the genuinely new Qi's. 

In this way we are guaranteed to find all the perfect forms. 
If we check for isometry with previously found forms the al- 
gorithm will at a certain point terminate, its output being the 
list of all perfect forms in a given dimension. The extreme 
rays of an n-dimensional polyhedral cone are the half-lines at 
which at least n~ 1 inequalities are binding (n — d{d + 1) /2 
here). The bottleneck of the algorithm is finding all the ex- 
treme rays Ri of a given lattice Q (or more rigorously of the 
Voronoi domain V(Q)), which, since the number of minimal 
vectors can be quite large (as much as exponential in d) can 
be a complicated linear programming problem. The generic 
version of this problem is known as a polyhedral representa- 
tion conversion problem in polyhedral computation commu- 
nity and its complexity is currently unknown"*' . All the 
forms generated from a given form Q are called neighbors 
of Q and the graph consisting of perfect forms linked to their 
neighbors is called the Voronoi graph of perfect forms in a 
given dimension d. Importantly, the graph is connected and 
starting from any vertex one can at least in principle reach any 
other vertex of the graph''*'''*''' '. 




FIG. 2. Example of lattice reduction for a square lattice: random 
initial basis {left) where basis vectors have large norms. After lattice 
reduction {right) one gets "short" basis vectors. 



Thus generated lattices might (and often do) have generat- 
ing forms with rather large norms of basis vectors. For ex- 
ample while we know there is just a single perfect form in 
d = 2. However a plain random walk would generate forms 
with entries growing as a function of step number To remedy 
this problem we use the fact that for a given lattice its basis 
can be transformed to an equivalent basis but with reduced 
basis vector norms. Figure 2 illustrates this idea for square 
lattice. The exact transformation which reduces the norms to 
the smallest possible value is expensive and we use a inexact 
one known as LLL-reduction after the names of the authors'*^ 
to produce equivalent representations of lattices with rather 
short basis vectors. Technically we apply the LLL-reduction 
on every newly generated form: this extra step allows us to 
generate forms with relatively small entries. Coming back to 
d = 2 case we find just 3 distinct forms (all of which are 
isometric). It is worth pointing out that the probability of gen- 
erating isometric forms becomes much less relevant for higher 
dimensions and completely irrelevant for d > 13. The LLL 



reduction is also a subset of isometry testing and actually re- 
moves the most trivial isometrics. 

In order to focus on higher dimensions we propose to ran- 
domize Voronoi's algorithm, namely to introduce a random- 
ized subroutine to find an extreme ray Ri. In this way we 
do not have to find all the extreme rays but just pick one and 
move in that direction. 

We do the following: we slice the cone with a plane, in this 
way the extreme rays become vertices of a polyhedron. We 
then define a random linear cost function 

d 

f{Q') - ^-oQ\, (26) 

where the Aij are gaussian random variables and we 
solve the corresponding linear programming problem 
maxQ/g-p(Q) f{Q')- Linear functions are necessarily maxi- 
mized at the vertices of the polyhedral region and therefore 
in this way we select randomly an extreme ray, which gives a 
neighbor of Q. The gaussian distribution of the Aij induces a 
distribution on the frequency each neighbor is visited which is 
far from uniform (a vertex is visited more often if, in the poly- 
hedron it is surrounded by facets with relatively large surface). 
We will discuss later our attempts to make more uniform this 
distribution. 

We have now defined the random generation of a new 
neighbor of Q so in order to define a random walk we need 
to define the rules for accepting or rejecting said moves. 

IV. MONTE-CARLO PROCEDURE AND THE VORONOI 
GRAPH 

It is clear that if we are only interested in the structure of 
the Voronoi graph we should run a random walk as unbiased 
as we can. Of course the most naturally unbiased algorithm 
would ideally generate any neighbor with equal probability. 
However this would be equivalent to finding all the neighbors 
for every perfect lattice; this problem can be incredibly diffi- 
cult and it has been solved only for d < 8 , with a large use 
of computer resources, so we do not attempt to solve it here. 

A. A warm-up: simple cases d < 7 

As a warm up we study very low dimensions: for d < 7 the 
problem of enumeration of perfect lattices is relatively simple 
due to small number of non-isometric perfect forms Af: 



Dimension 


1 


2 


3 


4 


5 


6 


7 




1 


1 


1 


2 


3 


7 


33 



The problem is completely trivial for d < 3 since there is a 
single perfect lattice (up to isometrics). For d ~ 4,5 enu- 
meration is trivial: our code finds the other forms on the first 
steps. Less trivial cases are d = 6 and d = 7 with 7 and 33 
perfect forms respectively. It takes about one thousand steps 
to find all 7 forms in d = 6. In d = 7 we recover 32 forms 
after 10^ steps. 



FIG. 3. Left The Voronoi graph in d = 6; vertex 1 is £5, vertex 3 is 
D(,, vertex 7 is, Aj. Right The Voronoi graph in d = 7: there are just 
33 perfect forms. The central point is Ej: it is connected to all the 
other vertices but A-j, which is the rightmost vertex of the graph. 



B. Properties of the d = 8 and d = 9 Voronoi graphs 

We compare the random walk on the exact Voronoi graph 
as found in with the numerical results of the previously de- 
scribed randomized Voronoi algorithm. 

The Voronoi graph for d = 8 is quite an interesting object 
if seen through the lens of statistical mechanics of random 
graphs. We unveil here only a small set of observations. The 
number of vertices is the number of perfect forms, namely 
10916, and we put an edge whenever two forms are Voronoi 
neighbors. The most connected form is the densest packing 
i?g, which has 10913 neighbors, and it is interesting to notice 
that the distribution of the connectivity of the graph follows 
quite closely a power law decay (a so-called scale-free net- 
work) for c < 2 10'^. Over this 3 -orders of magnitude range 
we can fit the connectivity distribution by the law 

p{c) oc c-(2-7±o.i)^ (27) 

which defines a critical exponent. We will see that this is also 
the case in d = 9. 

It follows from the large connectivity of E% that an unbiased 
random walk on this graph would visit E% a large number of 
times. By running a completely unbiased random walk on the 
exact Voronoi graph in 8 dimensions we find that £'§ should 
be visited about 1.6% of the times (this has to be compared 
with an average of 1/10916 ~ 0.01%). In our algorithm we 
see however that this number is much larger: is visited 
around 80% of the time. This means that our algorithm is bi- 
ased towards lattices with higher connectivity even more than 
an unbiased random walk is. This has to do with the large sur- 
face occupied by facets of the Ryshkov polyhedron enclosed 
by rays generating E^. 

This is a common feature in any dimension: the densest lat- 
tices are reached quite fast by our randomized algorithm even 
in absence of any a priori bias towards them. The balance be- 
tween the increase in the attractivity of the best packers and 
the increase in the size of the graph allows one to stumble 
upon the densest lattice up to li = 12 with a few hundred trials 
without having to bias the random walk towards the densest 
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FIG. 4. Top: The distribution of the connectivity of the d = 8 
Voronoi graph, exact results and Bottom: the same distribution sam- 
pled with the randomized Voronoi algorithm. AN{c) is the number 
of perfect lattices which have connectivity between c — 50 and c. 
The power-law fit is described in equation (27). In general, an un- 
derestimation by the random walk of the connectivity of the nodes is 
observed but a power law fit still works well, and the power law is 
compatible with the exact result (see text). 



lattices. Moreover, as a typical scale-free network, the diam- 
eter of the Voronoi graphs will be quite small, scaling as the 
logarithm of the number of vertices divided by the logarithm 
of the average connectivity. 

We now discuss the results of our randomized algorithm in 
d — We find, as said, that 80% of the times is spent on 
_Eg. The remaining 20% of the time is divided among the re- 
maining lattices. Every time a lattice is visited an isometry 
test is run against the previously visited lattices. If it is new, it 
is added to the list; in any case a link between the two lattices 
is added to the list of edges in the graph. In this way, in 10^ 
runs we generate around 3 • lO'^ non-isometric perfect lattices 
(out of 10916). This might be taken as a measure of the im- 
portance of isometry as well as of the dominance of iJg in 8 
dimensions. 

In cZ = 9 we run the randomized Voronoi algorithm for 10^ 
steps and we generate around 6 • 10^ non-isometric perfect 
forms. We recall that in d = 9 the Voronoi graph is conjec- 
tured to be made of around 2-10^ inequivalent perfect forms. 
We hence find in this case that the importance of isometry 
is much reduced. We will see that in higher dimensions the 
isometry test becomes irrelevant as randomly generated forms 
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turn out to be almost always non-isometric. 




ln(c) 



FIG. 5. The distribution of the connectivity of the d = 9 Voronoi 
graph estimated by the random walk. AN{c) is the number of per- 
fect lattices which have connectivity between c — 50 and c. The 
power-law fit is described in equation (28) 



By looking at the distribution of the local connectivity we 
see that also in this case a power-law distribution is the best fit 
over 3 orders of magnitude: 



p(c) cx c 



-2.5±0.1 



(28) 



We also observe the same slight overestimate of the fraction 
of low connectivity graphs we saw in d = 8. This is due (as in 
other dimensions) to the fact that the in order to assign a con- 
nectivity c to a graph the random walk has to visit said graph at 
least c times. There's no proved estimate of number of perfect 
lattices (size of the Voronoi graph) as a function of dimension. 
The sequence looks as 1, 1, 1, 3, 7, 33, 10916, ~ 2 • 10^ . . . 
and suggests a superexponential growth, for example like 
e . Consequently the number of steps required for an ac- 
curate estimation of connectivity grows rapidly. This means 
that for dimensions higher than 9 a different strategy has to be 
used. 

However, after observing the similarity between the two ex- 
ponents for the connectivity distribution and checking our ran- 
dom walk results against the exact results in d = 8 it is nothing 
but tempting to conjecture that the Voronoi graph is a scale- 
free random network in any dimension and that the exponent 
of the distribution of the connectivity is around 2.6. 

One can also plot (see fig. 6 and 7) the joint distribution 
of kissing number and energy observing how the best pack- 
ers have largest kissing number and they are both rare events 
with respect to the typical distribution. This phenomenon is 
constant across all dimensions. 



V. BIASING THE RANDOM WALK WITH A 
TEMPERATURE 

Following a common trick in statistical mechanics we intro- 
duce a temperature /3 as a Lagrange multiplier for the packing 
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FIG. 6. Top Kissing number vs energy {d — 8), generated set. Bottom 
Kissing number vs energy (d — 8), exact data. The insets show same 
plots with kissing numbers Z < 110. In both cases the best packer 
and kisser is alone in the upper left of the figures. 
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FIG. 7. Kissing number vs energy {d — 9). The inset shows detailed 
plot for kissing numbers Z < 140. 



fraction. We therefore would like to define a statistical ensem- 
ble described by the partition function: 



(Q) 



(29) 



e(Q) = -- log 0(g) 

where Q is a perfect lattice in d dimensions and is the 
measure induced on the space of perfect lattices by the solu- 



8 



tion of the linear program (26),^^ namely, ^J,{Q) is the fraction 
of times the lattice Q is visited when the random walk de- 
scribed in the previous section is run. We also defined energy 
of a packing e{Q) so in (10) that it is a quantity of order 1 
for the best packings which have packing fraction decreasing 
exponentially in dimension. Quite conveniently the best pack- 
ings translate into packings with lowest energy i.e. "ground 
states". The normalization for the temperature is due to the 
expectation that for the densest lattices log(0) ^ d (as both 
upper and lower bounds predicts) and we need the exponent to 
be order of the number of degrees of freedom, namely ^ (P. 

By lowering the temperature we expect to explore the re- 
gions of the Voronoi graph in which lattices are denser. 



VI. RESULTS 

Below we present the numerical results generated by ran- 
dom walks described above and their interpretation. 



A. Aims 

The generation procedure is inherently stochastic and we 
do not aim at generating complete sets of perfect lattices in a 
given dimension. As we already mentioned we have discov- 
ered 32 and approximately 3 • lO'^ forms after ^ 10^ runs in 
d = 7 and d — 8 respectively. The number of discovered 
forms in d = 8 increases with extra runs, although a complete 
enumeration would require a huge number of runs. 

Such a huge number of perfect lattices suggests a statistical 
approach so that properties of typical or even dense lattices 
can be extracted from a subset of the complete set. Thus our 
goal is rather to generate sufficiently large, representative sets 
of perfect forms in a given dimension which would allow us to 
understand typical properties of perfect lattices and spot any 
universal pattern behind. 

The fact that we are dealing with relatively large sets of 
forms together with the stochastic nature of the generating 
procedure allows to introduce empirical distributions of var- 
ious characteristics of lattices. We are going to focus mainly 
on two quantities: energy which was defined above and kiss- 
ing number. Both quantities are of interest with respect to 
the best packings. We will analyze their statistical properties, 
in particular their distributions and moments on the ensemble 
generated by the random walk. 

We have generated random walks (both simple and biased) 
in dimensions from 8 to 19. Complexity of computation grad- 
ually increases with dimension as does typical running time to 
generate sufficiently representative set of lattices. Runnning 
times vary from about an hour ind = 8,9to5 — 7 days 
in d = 19 to generate 5 • 10'* lattices. Higher dimensions, 
i.e. d > 20 are accessible, the difficulties encountered being 
rather of technical than conceptual nature. 



B. Random walk at infinite temperature 

We have first performed runs in different dimensions at in- 
finite temperature which correspond to plain random walks: 
departing from initial lattice one computes a random neigh- 
bour and hops there. It is natural to think that this way one 
generates typical perfect lattices . The walk terminates af- 
ter a finite number of steps N have been made. The averages 
(...) are simple summations normalised by N . 

Typically Ad was used as a starting point of a random walk 
for d < 12 and was used for d < 15—16 since the energy 
of Ad becomes too high. In even higher dimensions (d > 16) 
the energy of Dd itself becomes too high for Dd to be a good 
starting point and we used different initial lattices with better 
packing fractions which we generated by chain runs, that is 
first running a random walk starting at Dd and then picking 
a suitably dense lattice as a starting point for a new random 
walk. 

As already mentioned above our randomised code is biased 
towards denser lattices and it doesn't sample all lattices uni- 
formly like a complete enumeration would do (this effect is on 
top of the bias given by the larger connectivity of the densest 
lattices). It is instructive to compare our results to exact data. 
Unfortunately the latter are only known for d < 9 and there 
are too few perfect lattices for our approach to be benefitial for 
d < 8. So we start by comparing energy and kissing number 
distributions as sampled by our code and their exact values 
in d = 8. We see a reasonable agreement between the ex- 
act data and the ones generated by the randomised algorithm. 
This allows us to assume that data generated by randomised 
Voronoi's algorithm are representative and unbiased, and use 
data generated in higher dimensions where no exact data are 
available. The discrepancies present can be attributed to fluc- 
tuations associated to stochastic nature of our algorithm. This 
is especially clear for the kissing number which is integer by 
definition. 

A rough measure of representativity of a sample generated 
by a random walk is whether it visits "dense" lattices with 
high kissing numbers, or even better - the densest (known) 
lattice in that dimension. For low dimensions, d < 13, just 
Nd = 10* runs were enough to satisfy this requirement. Start- 
ing with d = 13 one has to make more runs (Although in 
d = 13 a random walk of 10"* steps comes quite close to the 
best packing: e = 0.28 and Cbest ~ 0.27). The required num- 
ber of steps Nd is growing fast: iVia ^ 10^, while Nu > 10^. 
The situation quickly deteriorates in higher dimensions: while 
in d = 8 random walk is hitting Es around 80% of the time, 
the number drops down to < 1% of hits for Aio - the best 
known lattice pakcing in d = 10 - and goes further down for 
higher d. Table I gives frequencies for a random walk to visit 
the best packers in d = 8 — 12. The data seem to suggest a 
faster than exponential decay, a simple fit giving ^ e~^ ° ^ 

The table II gives a summary on average energies, their 
standard deviations (Te, best found, worst found and best 
known lattice for d = 8 — 19 {N is number of steps in random 
walk): Standard deviation clearly decreases with dimensions; 
the increase for d — 17 — 19 indicates that more runs are 
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FIG. 8. Comparison of exact and empirical distributions generated 
by randomised Voronoi's algorithm. Top Distribution of energies e 
from randomised Voronoi's algorithm with isometry testing (blue) 
and exact distribution (red) for d = 8. Bottom Same comparison 
of distributions of kissing numbers from randomised Voronoi's al- 
gorithm with isometry testing (blue) and exact distribution (red) for 
d = 8. 



Dimension 


8 


9 


10 


11 


12 


Frequency 


0.835 


0.341 


0.096 


0.0156 


0.00191 



TABLE I. Frequencies with which a best known packing is visited 
by a random walk as a function of dimension. 

required to get a representative set of lattices. Indeed compar- 
ing behaviour of the deviation with number of runs for d = 17 
one sees the decrease as number of runs increases (the same 
behaviour is present in d = 18, 19): The decrease of standard 
deviation suggests that distribution of energies Vd{e) is con- 
centrating around mean value and becomes peaked around its 
mean value for large d and for d = cx): 

Vd^oo{e)-5{e-{e)d^^). (30) 

Fig. 10 shows behaviour of average energy (no checks for 
isometry) with dimension. Large deviations in low dimen- 
sions up to d < 12, represented by errorbars on the figure, 
are related to the fact that the distribution of energies in these 
dimensions is highly irregular if no check for isometry is per- 
formed during the random walk (see Fig. 9, case of d = 8 for 
an illustration). 

An important issue is equivalence/isometry of generated 



Dim. 


(e) 




Best found 


Worst found 


Best known 


8 


0.180572 


0.021502 


0.171465 


0.308792 


0.171465 


9 


0.23352 


0.01823 


0.21396 


0.34188 


0.21396 


10 


0.26828 


0.01422 


0.23857 


0.37285 


0.23857 


11 


0.29347 


0.01228 


0.25511 


0.40193 


0.25511 


12 


0.31505 


0.00796 


0.25055 


0.38024 


0.25041 


13 


0.33106 


0.00328 


0.27179 


0.40709 


0.27178 


14 


0.34277 


0.00236 


0.31862 


0.43265 


0.27386 


15 


0.35405 


0.00273 


0.33522 


0.45703 


0.27218 


16 


0.36507 


0.00197 


0.34235 


0.48031 


0.26370 


17 


0.37205 


0.00280 


0.33949 


0.50258 


0.27833 


18 


0.38322 


0.00235 


0.37805 


0.39238 


0.28489 


19 


0.39000 


0.00391 


0.37909 


0.40146 


0.28903 



TABLE II. Average energies of perfect lattices for d — 8 • ■ • 19. 
Sample sizes N are 10^ for d = 8 - 12, 10^ for d = 13 - 16, 
2 ■ 10^ for d = 17, 18 and 1.5 • 10^ for d = 19. The observed 
increase of standard deviation for d > 17 indicates that sample 
size was not big enough. Increasing the sample size decreases the 
deviation. 



iV 


10^ 


10^ 


2 • 10^ 


Deviation 


0.0064997 


0.0030565 


0.0028058 



TABLE III. Standard deviation of energy in d = 17 as a function of 
number of runs A'^. 



lattices. As we have discussed above a single lattice admits 
many equivalent representations in terms of quadratic forms. 
One might worry if random walk is generating many/few 
equivalent lattices. The above results were generated neglect- 
ing isometry partially: ony LLL-reduction was performed on 
newly generated forms. Based on d = 7, 8 results we know 
that isometry is definitely important in low dimensions. How- 
ever it is only relevant for low dimensions, our data suggest 
d < 13, where the number of perfect lattices is relatively small 
and random walks of moderate size contain many isometric 
copies of the same lattice. For higher dimensions, d > 13, 
where the number of perfect lattices is huge the chance of hit- 
ting an isometric lattice is vanishingly small except for the 
densest lattices which have a larger isometry family. This is 
illustrated by Fig 9 which compares probability distributions 
of energies for d = 8 and d ~ 12. It is worth stressing that 
this statement holds true only if one samples a relatively small 
subset of all perfect lattices. Once sample size is comparable 
to the size of the full set of perfect forms, isometry becomes 
important in any dimension. This fact can in principle be used 
to define a formal criterium whether one has generated a re- 
spresentative sample. Including isometry test in generation 
procedure is easy: every newly generated form is checked for 
isometry against all previously genrated forms . 

The effect of isometry on energy average ( e) is to in- 
crease values for low dimensions, which are dominated by 
dense lattices if no isometry cheks are performed. The higher- 
dimensional data, d > 12 are left intact since isometry be- 
comes completely irrelevant. We reproduce the table IV and 
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FIG. 9. Top. Distribution of energies e with (blue) and without 
(red) isometry test in d = 8. Isometry is very important and the two 
distributions are completely different: without isometry check distri- 
bution concentrates around the energy of _Eg). Bottom. Distribution 
of energies e with (blue) and without (red) isometry test in d = 12. 
Isometry is no longer important and the distribtuions are almost the 
same, except low energy tail, where one still sees small spikes. 



d 


(e) 


Std(e) 




Std,(e) 


(e)ex 


Stde,(e) 


8 


0.180571 


0.021502 


0.258296 


0.0050364 


0.258845 


0.003593 


9 


0.233521 


0.018231 


0.266341 


0.005073 


0.259662* 


0.006006* 


10 


0.268281 


0.014227 


0.281615 


0.005484 






11 


0.293471 


0.012288 


0.299142 


0.005262 






12 


0.31506 


0.007967 











TABLE IV. Comparison of energy averages { e) without and with 
isometry test. Additionaly exact values of average and standard de- 
viation are given for d = 8. * Values for d = 9 are extracted from 
partial enumeration . 



the ( e) curves for data with isometry checks. We see the same 
trend of decreasing standard deviation with increase of dimen- 
sion as in the case of no isometry testing. In what follows we 
are using samples with checks for isometry for d < 12 and 
with no isometry checks for d > 11. 

In what follows we use mixed set of data: samples with 
isometry checks for d < Yd and samples with no isometry 
testing applied for d > 12. We do so to remove features spe- 
cific to low dimensions d < 13 and reveal the generic features 
common with dimensions d > 12. 



FIG. 10. Average energy (e) = —^(logip) of a random walk as a 
function of dimension (d — 8—19). Errorbars correspond to standard 
deviation of energy. The smooth curve is a guide to the eye. 




FIG. 11. Average energy (e) = — ^(log^) of a random walk as a 
function of dimension (d = 8 — 19) (red crosses) compared to en- 
ergy of Ad and Dd lattices (see Eq.(31)). The dashed curves are the 
leasing asymptotics of Minkowksy (top), Torquato-Stillinger (mid- 
dle) and Kabatiansky-Levenstein (bottom) bounds. The Minkowksy 
and Torquato-Stillinger are upper bounds while the Kabatiansky- 
Levenstein bound is a lower bound on the energy of the best packing. 
The yellow and blue continuous lines are Ad and Dd, the green and 
violet lines are the two fits (37), (38). 



Let us concentrate on two possible scenarios, the simplest 
cases where to locate our typical lattices. On one hand we 
can for example look at the energies of Ad, Dd families of 
lattices : 



e{Ad) 



e{Dd) 



2 log - + 
log(d)/2 - 

log(d)/2 - 



log(l + d) , 1 



2d 



0{1) 
0{l) 



(31) 
(32) 



-)log2+-logr(l + -(p3) 



Both Ad and Dd have asymptotically equal energies for d 
oo: ^ log d/2 which means sub-exponential packing fraction. 

Minkosky's and Kabatiansky-Levenstein bounds tell us that 
there are lattices with only exponentially small packing frac- 
tion. Asymptotically in large dimensions, upper and lower 
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bounds give: 

eA./ = log(2) + 0(log(d)/d) (34) 
e^i = 0.413... (35) 

and it is worth remembering the Torquato-Stillinger conjec- 
tured bound which should replace Minkoswky's under appro- 
priate hypothesis on high-dimensional lattices'"'^: 

cts = 0.539 + 0(log(d) /d) . (36) 

Random walks in high dimensions are sampling lattices 
with energy close to its mean value (e). We try two fits for 
this function of d, one with the leading order term constant, 
hypothesizing a "best packer" behavior for typical lattices in 
high dimensions and the other with leading log(d)^' . For the 
first we obtain 

(e) = (0.58 ± 0.04) - ^°4^(0.9 ± 1.0) - (0.8 ± 0.6)^-^. 
d 

(37) 

The constant term is suggestively close to the Torquato- 
Stillinger bound and, within the associated error, it is below 
the Minkowsky bound log(2) = 0.69. However, an equally 
good fit can be obtained by assuming that the leading term is 
growing logarithmically 

(e) ^ (0.066±0.04)log(rf) + (0.27±0.04)-(1.4±0.2)d-i 

(38) 

although the coefficient of the logarithm is well below the 
value 0.5 of the Ad and Dd families (typical lattices are much 
denser than these examples). Both fits are equally good, as 
can be seen from Fig. 10, the resolution of the two can only 
occur for d ^ 40. 

The main effect of isometry on distribution of energies 
7^(e) is to supress low energy spikes (see Fig. 9) associated 
with dense lattices which are relatively often visited in these 
dimensions by a random walk, and shift the weight to the uni- 
versal bell-like feature which dominates the distribution Vdie) 
in high dimensions. As of the distribution of kissing numbers 
Z switching on the isometry testing kills the large-Z tail of 
the distribution and concentrates the weight around small val- 
ues of Z of order d{d + 1) (recall that this is the lower bound 
on kissing number for perfect lattices). These facts indicate 
that in high dimensions typical perfect lattices have relatively 
high energy (but still lower than Ad and Dd) and small kiss- 
ing numbers, of order d{d + 1). If we define rescaled vari- 
able X = {e — {e)d) /cTe we expect the probability distribution 
functions of x to collapse on some master curve with mild 
dependence on d: 

Vd{x) oc Vd (^^) ■ 

Indeed after rescaling a master curve is emerging as shown 
on Fig. 13 though the collapse is not perfect: case d = 12 
is special with quite different shape as compared to other 
dimensions as highlighted on Fig. 13. All the distributions 
are skewed to the left, i.e. towards denser lattices, although 
this is hard to spot on Fig. 13 while this is clearly so for 
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FIG. 12. Probability distributions Vd of energy e for d = 8—10, 12— 
19 - color goes from red (d — 8) to violet (d — 19). As dimension 
increases averages increase and peaks shift to the right. 



d = 12. These features become more pronounced if one 
studies gd{x) = — log Vd{x) showed on Fig. 14: the generic 
skeweness to the left (towards the denser lattices) becomes 
clear. For all dimensions studied except c? = 12 the central 
part of gd{x) can be well fited with a Gaussian 

-\ogVd{x^Q)^O.^h+^, 

i.o 

the value of the coefficient of x"^ being slightly larger than (but 
still consistent with) 1/2 reflects the skewness of the distribu- 
tion. The skewness only appears for larger values of x which 
are noisy because we do not have enough statistics to probe 
them accurately. 




FIG. 13. Le/f Probability distributions 75(2;) ford = 8-10, 13-17- 
color goes from red for d = 8 to magenta for d = 19. We have used 
exact distribution for d = 8 for convenience and skipped d = 12. 
Right Comparison of distributions Vd{x) for d = 10, 12, 13. The 
case d = 12 is very disctinct from neighbouring dimensions. 

We now study the statistics of kissing number. For a typical 
perfect lattice the kissing number is of order d^, i.e. like for 
Ad or Dd, and of the same order of magnitude as the lower 
bound d{d + 1). To highlight this point we normalized {z) 
by d{d +1), the minimal possible kissing number which gave 
a curve shown on Fig. 16. Thus a typical perfect lattice is 
similar to Ad or Dd in kissing numbers but has a lower en- 
ergy/higher packing fraction. As we see from Figs. 15 and 
1 6 kissing number fluctuates much stronger than energy and 
the only conclusion we can make from the plots is that the 
distributions concentrate around their means just like it hap- 
pens with energy. Combining this observation together with 
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FIG. 14. Gaussian fit to the central part of the probability ditribution 
7'd(x)ford = 8-ll,13-19. 
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FIG. 16. Average kissing number normalised by ci(d+l) ford = 8 — 
19. Errorbars correspond to first and third quartiles (These are zero 
for d — 18, 19). Despite strong fluctuations the value of normalised 
kissing numbers is of order 1. 



behavior of average energy we see that in high dimensions the 
Voronoi graph is dominated by lattices which have properties 
similar to Ad and Dd- 
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FIG. 15. Average kissing number for d = 8 — 10, 12 — 19. The red 
curve is the best known kissing numbers in corresponding dimen- 
sions. 



C. Random walk with ^ > 

As dimension is increased beyond d ^ 13 we are no longer 
able to recover the densest known lattice packing with a plain 
random walk, at least for the number of steps we have tried 
(from a few hundred thousands to a few millions, depending 
on dimension). Given a fast growth of the number of peifect 
forms with dimension, one would likely have to sample ran- 
dom walks of size comparable to the number of perfect forms 
to see the densest lattices, something that is out of reach al- 
ready for moderate dimensions c? ~ 13 — 14. 

We therefore introduced a procedure which biases the walk 
towards denser lattices. We employed standard Metropolis- 
like rule with fictitious temperature (3 described above in 
Sec. V which favours denser lattices. Namely, we generate 



a neighbor Q' of the lattice Q and compute its packing frac- 
tion 4){Q') and from this its energy e{Q'). If e{Q') < e{Q) 
we accept the move and if e(Q') > e{Q) we accept the move 
only with probability exp{~f3{e{Q') — e{Q))). 

This allowed us to recover consistently the densest (known) 
lattice packings up to d = 17 and to get very close to the best 
known lattices in d = 18, 19, where we start seeing some 
complex landscape behavior. We managed to get the best 
known pakcing in these dimensions too but in a much less 
consistent fashion. 

Again we are looking at distributions and moments - aver- 
age and standard deviation - of energy and kissing number 
We saw for plain random walk which coiTesponds to (3 — 
that E{d) = ( e) is a smooth curve as a function of dimen- 
sion. As the temperature is lowered Ep{d) curves become 
more singular reflecting the peculiarities of any given dimen- 
sion: it is well known that the nature of dense sphere packings 
varies greatly as a function of dimension - one of the factors 
that makes the problem of sphere packing so complicated. Up 
to d = 11 changing the temperature immediately affects the 
range of energies probed by the random walk: the lower the 
temperature the lower the energy and E{/3) — ( e)^ is essen- 
tially an exponentially decaying function of /3. Starting from 
d = 12 and up the pattern of E{f3) changes qualitatively: a 
plateau emerges at small (3 where the probed energy is almost 
insensitive to variations of temperature and is roughly equal 
to energy of /3 ~ random walk. As inverse temperature 
/3 is increased there is a crossover to lower value of energy. 
The value E{f3) for large /3 is approximately equal to the 
ground state energy, again almost insensitive to variation of 
/3. Furthermore, sufficiently close to the crossover we observe 
strong run to run fluctuations of values of ( e)^, a phenomenon 
which is reminiscent of a glassy free energy landscape . Such 
behavior suggests a phase transition as a function of /?, as 
d — > oo: as the temperature is lowered one leaves a universal 
phase dominated by typical perfect lattices and enters a phase 
where lattices with low energies dominate the biased random 
walks. To test this assumption we define f3c{d) as a solution to 
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Edipc)) = E,{d) = {EM + Ed{^))/2. As usual Ed{^) 
should read as Ed{(3i) for some sufficiently large The 
crossover width is defined as /3^{d) — P>{d) where 



Ed{0) - Edi^) 



E<{d) = Ed{p<) = Ed{^) + ^Arf = -^Ed{0) - \Ed{^) 
Ey{d) = EdiP>) = Edic^) + ^Ad = ^Ed{0) - ^Edic^) 



FIG. 17. Top Average energy (e) — — |(log<^) of a biased random 
walk as a function of dimension d = 8 — 17: inverse temperature /3 
goes from (red) to 5 (violet); red crosses are the best known lattice 
packings. As the temperature is decreased, details of the scenarios in 
finite dimensions become relevant. 




FIG. 18. Average energy (e) = — ^ (log (p) of a biased random walk 
as a function of temperature for dimensions d = 8 — 1 1 (color goes 
from red to green); dashed lines are the best known energies in cor- 
responding dimensions. 



The choice of factors 1/4 and 3/4 is not important and they 
can be replaced by other number. If there is indeed a phase 
transition then W = {Py— li<) / Pc should converge to a con- 
stant value as d — > oo. Fig. 20 shows dependence of W on 
dimension. One observes indeed a tendency to convergence 
to a constant value of 0(1) (although with noticeable oscilla- 
tions around it). We attribute the increase for d > 17 to the 
glassy nature of the energy landscape of perfect lattices: these 
are exactly the dimension where the simple Monte-Carlo ap- 
proach starts experiencing problems finding the best packer. 
The ci = 18 is intermediate between d < 18 and d = 19. 
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FIG. 20. = (S> - S< ) /Ec as function of dimension d = 8 - 18. 
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FIG. 19. Average energy (e) = — ^ (log ip) of a biased random walk 
as a function of temperature for dimensions d = 12 — 17. 



The situation seems to change qualitatively in d = 19: for 
mildly low temperature one has to increase drastically run- 
ning time (as compared to d < 18) in order to reach the best 
known packings. For very low temperatures, /3 ^ 3 — 5 for 
d ~ 18, 19, the Monte-Carlo routine gets stuck around some 
relatively dense lattices and is never able to recover the dens- 
est lattice, or even approach it within the accuracy achieved in 
smaller dimensions. Typical energies reached by Monte-Carlo 
are of order e ~ 0.35 — 0.36 for /? < 5. This is to be compared 
to the ground state e = 0.29 coiTesponding to lattice Aig. It 
is then crucial to study higher dimensions in order to under- 
stand whether this behavior is a peculiarity of d = 19 or it is 
a generic trend establishing in high dimensions. However we 
are unfortunately currently unable to investigate dimensions 
higher than 19 but we hope to be able to do so in the future. 
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VII. DIAMETER OF THE VORONOI GRAPH 

An interesting question is the number of perfect forms as a 
function of dimension d. The exact numbers for d < 9 and the 
estimate in c? = 9 suggest very steep, perhaps superexponen- 
tial law which would make the full enumeration impossible 
beyond d ^ 11. We conjecture that the number of perfect lat- 
tices should grow as JVd ^ cxp{A d^) for an appropriate con- 
stant A for large d. This conjecture is natural in the framework 
of statistical mechanics as the number of degrees of freedom 
is 0((i^) and so should be the "entropy" of the system. 

Looking at the distribution of the coefficients we can more- 
over conjecture that the Voronoi graph is a scale-free random 
graph, at least for a range of connectivities and for large d. 
For scale-free networks an estimate of number of vertices as a 
function of connectivities c of the vertices of the graph is - 



(log c) 



(39) 



Here Diam(CJd) is diameter of the graph: the longest among 
the shortest paths between any pair vertices. 

We have estimated the diameter of the Voronoi graph Qd 
using the information on the graph provided by the random 
walk. This contains partial information and serves just as an 
order of magnitude consideration so we must consider the de- 
pendence on the size of the sample. This computation be- 
comes increasingly harder with growing d and we have re- 
stricted the study to d < 11. 

If the distribution of the connectivity is indeed scale free 
with fixed exponent 2.6, we find that 



(lege) = — — = 0.62, 
\ 6 / 2.6-1 



(40) 



We find a reasonable agreement with numerical estimates of 
(lege): 1.274,0.954,0.771,0.7 for d = 8, 9, 10, 11 respec- 
tively. The excess of values of (logc)^ with respect to con- 
jectured value 0.62 is due to the fact that we sample many 
well connected, dense lattices while not visiting many lattices 
with low connectivity. Therefore the logarithm of the size of 
the graph and the diameter should be proportional as 



logA/'<j~0.62Diam(ad). 



(41) 



We can then test if our hypotheses on the connectivity, the 
number of forms and size of the graph fit well together 
We find graph diameters 3, 6, 13, 32 and 131 for d ~ 
7, 8, 9, 10, 11 respectively. Remark that the exact diameter is 
3 and 4 in d = 7 and 8 respectively. The growth is clearly 
faster than linear as shown on Fig 21 and is consistent with 
the hypothesis of scale-free Voronoi graph. Quadratic fit for 
Diam{Qd) based on data for d = 7 — 10 reads as: 

DmrniGd) = 217.6 - 58.6 d + 4 d^ 

However with the actual data we cannot find the precise 
scaling. Although exponential fit looks more accurate than 
quadratic on Fig 21 we know that there are many forms in 



d = 11 which were not visited by a random walk. Their addi- 
tion to the graph would reduce the diameter and perhaps smear 
the seemingly exponential growth. More data are required to 
resolve this issue and we leave the resolution of this problem 
for future work. 




FIG. 21. Red crosses: estimate of diameter of the Voronoi graph 
as function of dimension. Blue and brown curves are quadratic and 
exponential fits respectively provided here as guides for the eye. 



VIII. TRYING TO UNIFORMIZE THE CHOICE OF 
NEIGHBOR 

As we have already mentioned above, the randomization 
of Voronoi's algorithm is not unique: different cost func- 
tions (26) produce slightly different results. We have consid- 
ered a number of functions, targeting uniformization, i.e. try- 
ing to make sampling of rays/neighbors more uniform, more 
like it is for full enumeration. In all cases we observed a bias 
towards denser forms with higher kissing numbers, which we 
try to reduce. In particular we constructed a "uniformized" 
cost function as shown on Fig. 22 (recall that we have a n- 
dimensional polyhedron, n = d{d + l)/2 here, defined by a 
set of inequalities, the number of inequalities J\f > n). This 
construction is inspked by the remark that purely random cost 
function generates rays weighted with areas of facets adjacent 
to that ray, and it also favors forms that have higher connec- 
tivity, i.e. number of neighbors. This is an advantage if one is 
interested in denser forms. However if one is studying prop- 
erties of the Voronoi graph it might be preferable to make the 
outcome of neighbor generation more uniform. 

The above construction tries to give facets a more uniform 
weights. Comparison of numerical results for random and uni- 



Input: Voronoi domain V((5) 

Pick an inequality at random 

Saturate the inequality, i.e. replace it with equality 

Make a random Gaussian cost function / as before 

Solve linear program to get an extreme ray 
Output: Random extreme ray R 

FIG. 22. Algorithm for uniformized random extreme ray generation. 
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form cost functions are presented on Fig. 23 which shows dis- 
tributions of kissing number and energy in d = 8. There's 
no significant difference of distributions between the random 
and unformised cost function. However the uniformized cost 




0.22 0.24 0.26 0.28 0.30 

e 
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FIG. 23. Top. Distribution of energies Bottom. Distributions of kiss- 
ing numbers. Blue and red curves are generated by random walks in 
d = 8 with random and uniformized cost functions. 



function is advantageous over the random function if one is 
interested in the properties of the Voronoi graph: typically it 
yields more non-isometric forms than the pure random func- 
tion for equal number of runs. We have performed this com- 
parison for d = 8 — 12 and results are summarized in the table 
below (where Fraction column is a ratio Afu/J^r of number 
of forms found Mr and Mu with random and uniformized cost 
functions respectively): 



Dim. 


Steps 


Random 


Uniformized 


Fraction 


8 


210^ 


1793 


2955 


1.648 


8 


410^ 


2529 


3963 


1.567 


10 


10^ 


331065 


434317 


1.312 


11 


10^ 


744282 


825695 


1.109 



The difference between the two cost functions is decreasing 
rapidly as dimensionality is increased. We believe that these 
strategies are better suited for lower dimensions d < 12 where 
isometry is important. 



IX. EXTENSION TO PERIODIC SETS 

Before concluding let us describe a possible extension of 
our approach to lattices with many particles per unit cell 
which we refer to as periodic sets throughout this section. 
Such an extension is possible but has a number of limitations 
which make the problem more difficult than the Bravais lattice 
version. 

The generalization of the Voronoi algorithm to periodic sets 
was introduced by Schiirmann . An m-periodic set is defined 
by a quadratic form which describes how a unit cell is trans- 
lated in space and a set of m real vectors {translational part) 
that defines the positions of m particles inside the cell. It is 
then possible to extend the Voronoi theory presented in Sec. 11 
and introduce m-perfect and m-eutactic lattices; m-extreme 
lattices are defined as local maxima of packing fraction of m- 
periodic sets, just like in the Bravais case. There is as well an 
analogue of the Ryshkov polyhedron. 

It is at this point that a crucial difference appears which 
makes the problem more complicated than the lattice one. In 
general not all extreme lattices are m-perfect and 771-eutactic: 
there exist lattices which are extreme, but not perfect. An ex- 
ample is provided by fluid diamond packings where a fraction 
of spheres can be moved around freely without canghing the 
packing fraction. Furthermore the Voronoi graph no longer 
exists: the method only provides a local direction in which 
packing fraction is increasing. Potentially this allows to de- 
sign an algorithm that starts with a periodic set and end up 
at an m-perfect lattice . On the other hand, the extension to 
many particles in a unit cell highlights the importance of per- 
fect, strongly eutactic lattices since one can prove that they 
are extreme ", that is they are extreme among lattices with any 
number of particles per unit cell. 

These limitations are lifted if one fixes translational part and 
replaces real vectors in the definition of a peridic set by their 
rational approximations . Under this assumption, all the fea- 
tures of the Voronoi theory are recovered. Yet the complexity 
is increasing too: the computation of the shortest vectors of 
such periodic set is more involved. 

X. CONCLUSIONS AND FURTHER DIRECTIONS 

We have suggested a new approach to the lattice sphere 
packing problem based on randomization of the Voronoi al- 
gorithm. Previous works used complete enumeration that be- 
comes computationally unfeasible beyond d ~ 10 — 11 (see 
however '^'^'"^^). We have developed an implementation of 
our algorithm that allowed us to study dimensions from 8 to 
19 and we foresee its application for studying perfect lattices 
up to d = 40 at least (beyond that, technical problems with 
the implementation of the algorithm become conceptual prob- 
lems). 

We have studied statistical properties of the sets of perfect 
lattices generated by our algorithm, both typical and extreme 
values focusing on two quantities: energy, which we define as 
proportional to the logarithm of the packing fraction, and the 
kissing number. For all dimensions except d = 19 we were 



able to retrieve the best known packings starting from Ad or 
Dd lattices either using simple random walk for d < 12 or 
biasing the random walk with temperature for d > 12. In 
d — 19 we had to restart the walk many times in order to 
hit the best packer: random walk was always getting stuck 
in some higher-energy lattice, a phenomenon which is remi- 
niscent of a glassy free energy landscape. The change of the 
average energy with temperature suggests the existence of a 
sharp phase transition as d — > oo, although we cannot argu- 
ment on this topic more, due to the large dimension-dependent 
fluctuations as the energy is lowered. We do not exclude we 
will be able to say more on this topic in future work. 

We also found that the typical values tend to have much 
smoother behavior what allowed us to propose two possible 
scenarios for the large d behavior of the packing fraction of the 
typical perfect lattices: in one case we obtain en exponential 
decay of the packing fraction whose leading order improves 
upon Minkowsky's bound 

_ 2-(0-»4±"-06)^ (42) 

while in the second case we have a faster, factorial-like decay 

however with an unnaturally small exponent. The resolution 
of this conundrum would need investigation of lattices in di- 
mensions 40 and higher. 

Higher dimensions are also accessible and will require 
mostly technical rather than conceptual modifications in the 
code, at least for d < 40. Getting beyond d = 24 is quite 
important since in dimensions below 24 are dominated by the 
Leech lattice A24 and all the densest lattices in these cases are 
cross sections of A24. 

Other possible applications of our work include a test of 
the "decorrelation principle" in , by studying the two-particles 
correlation functions of typical perfect lattices, and a system- 
atic study of the perfect and eutactic lattices which are the 
true local minima of the energy for the purpose of unveiling 
a glassy structure of the energy landscape. Checking for eu- 
taxy is quite straightforward, after a set of perfect lattices has 
been generated, but we found that this requires a much larger 
statistics than that used in our paper since the rejection rate is 
quite large: as dimension of space is increased the fraction of 
(at least) eutactic lattices discovered by a plain random walk 
drops rapidly as illustrated in Table V. If one biases the walk 
with temperature the numbers increase, but they are still low 
and we have not tested whether the increase is due to different 
lattices or isometric copies of few lattices. Therefore we leave 
this for future work. 

Finally, randomization procedure we have introduced could 
also be applied to other optimization problems like lattice cov- 
ering problem , where one searches for the most economical 
way of covering a space with spheres of equal size. Another 
possible activity along the same direction is to adapt our ran- 
domization procedure to the algorithm generating all eutactic 
lattices in a given dimension . 

As we have indicated, finding extreme rays of the Voronoi 
domain V is a particular case of a general polyhedral repre- 
sentation conversion problem . This is an important problem 
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Dimension 


Fraction of eutactic lattices 


8 


0.997 


9 


0.830 


10 


0.738 


11 


0.479 


12 


0.134 


13 


5.11e-03 


14 


3.00e-04 


15 


1.30e-04 


16 


8.00e-05 


17 


6.00e-05 


18 


1.25e-05 


19 


2.00e-05 



TABLE V. Fraction of eutactic and strongly eutactic discovered by 
random walk for d = 7 — 19. 



in combinatorial optimization and computational geometry. 
Although efficient algorithms exist for certain classes of poly- 
hedra, its complexity in general is unknown • but all ex- 
isting algorithms, that perform the full conversion, are expo- 
nential in dimension of a polyhedron . In this wider context 
our randomization approach offers a possible workaround for 
optimization problems which require solution of the represen- 
tation conversion problem in order to find an optimum. 
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APPENDIX A. SOME TECHNICAL DETAILS 

The two main techincal ingredients of the Voronoi algo- 
rithm are generation of random extreme ray R of the Voronoi 
domain V{Q) and finding a neighbour Q' of a given lattice Q 
provided an extreme ray R. 

Computing a random extreme ray has the same complexity 
as generating the Voronoi domain V{Q) and solving a linear 
program. We need to know shortest vectors of Q in order to 
build V{Q). Computing shortest vectors of a lattice is expo- 
nentially hard problem in d. However decent algorthims exist 
allowing computation to caried out in reasonable time at least 
up to d ^ 40 ■ \ The other source of complexity is the size 
of linear program which is defined by kissing number of Q 
(and hence scales exponentially in d for dense packings) and 
is limited by ability of linear program (LP) solvers to cope 
with huge linear programs: size of LP becomes of order 10^° 
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Input: perfect form Q, extreme ray R 

while Q + uR^ SU and \(Q + u R) = \{Q) do 
itQ + uR^ Sio and A(Q + u i?) = A(Q) then 

«^ (; + m)/2 

else 

u) <— (ii, 2u) 

end if 
end while 

while Min(Q + I R) C Min(Q) do 

g^{u + l)/2 

it\{Q + gR) > A(Q)then 

I ^ g 

else 

u ^ min{(A(Q) - Q[v])/R[v]\v £ Mm{Q + 
gR),R[v] < 0}U{.g} 

end if 
end while 
Output: a I 

FIG. 24. Modified binary search for neighbour Q' of a lattice Q 
given an extreme ray R. 



for the densest known lattices in d > 40. Based on this obser- 
vations we expect our method to work up to d ^ 40, at least 
in theory. It is also worth pointing that it is straightforward to 
check if a given ray R is extreme . 

Finding a neighbour Q' ~ Q + aR with a G Q proved 
to be a harder problem computationally and it is this part 
of the problem that put limited our data by d < 20. The 
value of a is rational ■ , so that we can always choose Q' 
to be integral and all perfect lattices then have integral rep- 
resentation. We use modified binary search algorithm as de- 
fined by Schiirmann to compute neighbours of a lattice {S^q 
is set of all lattices) presented on Fig. 24. The idea be- 
hind this construction is very simple: the neighbour of Q is 
Q' — Q + a R with the smallest positive rational a such that 
A(Q) = A(Q + a R) and Min(Q + aR) % Min(Q) In the 
first part above upper and lower boundaries for a are defined. 
The second part is a modified binary search for value of a. 
The modification - an extra conditional in the assignment of u 
- is necessary to make the algorithm converge in finite number 
of steps to an exact rational value of a. 



APPENDIX B. RANDOM WALKS AND ISOMETRY CHECK 

We have used two different approaches to perform checks 
for isometry of lattices. In the first approach we split the data 
generation in two steps 

• Generate a random walk in space of lattices with no 
check for isometry. 

• Run isometry test on the trajectory of the random walk 
and generate an approximate Voronoi's graph. 

After the first step one obtains a full trajectory of a random 
walk as list of lattices. The second step generates the graph 
by eliminating isometric copies of lattices by glueing together 



Input: perfect Q, graph G = (V = 0, £ = 0) 
loop 

Random extreme ray R -i^ Q 
Neighbour Q' -(^ Q + aR 
for P e G do 
if P ~ Q' then 

E ^ EU{Q,P) 

else 

V ^VUQ' 
E^EU{Q',Q) 
end if 
end for 
end loop 
Output: Voronoi graph G 

FIG. 25. Algorithm that constructs an approximation to the Voronoi 
graph 



Input: Lattice A 
loop 

(*) A' ^ A 

Accept A' with some probability p 
Goto (*) 
end loop 
Output: Dense lattice A 

FIG. 26. Naive random wallc 



isometric elements of the list. This induces a relation of neigh- 
bourhood in the list and transforms the list into a graph. Sec- 
ond possibility is to perform isometry check and graph con- 
struction on the fly (P ~ Q denotes isometric equivalence, 
V and E are sets of vertices and edges of the graph G re- 
spectively) as shown on Fig. 25. Algorithm terminates after a 
predefined number of steps has been done. 

An algorithm to check whether two lattices are isometric 
was developed by W. Plesken and B. Souvignier in Ref. . We 
adapted the original code of B. Souvignier to perform isome- 
try testing. 



APPENDIX C. NAIVE RANDOM WALK 

It is worth discussing performance of a straightforward ap- 
proach one might be tempted to follow. The Voronoi con- 
struction is elaborate and requires computational effort. A pri- 
ori one might wonder if a simple lattice random walk/Monte- 
Carlo is preferrable (maybe in higher dimensions) ? The al- 
gorithm shown on Fig. 26 is extremely simple: one hopes to 
approach the best packer by small steps if the random walk is 
sufficiently biased towards denser lattices. When generating a 
move one has the option of eigther producing a new lattice A' 
which might or might not be an isometric copy of A. Accep- 
tance probability p could be 1 (random walk) or for example. 
Metropolis rule (a la Monte-Carlo). 

Unbiased random walk (infinite temperature in our lan- 
guage) with moves that generate non-isometric lattices 
A' gives an average packing fraction which is equal to 
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Minkovsky's bound^""''^ This is a rather strong result since 
Minkovsky's bound is non-constructive and constructing a lat- 
tice in a given dimension satisfying the bound is yet an open 
problem. However it is very hard to implement that type of 
updates in practice " and one has to rely on various approx- 
imations. In case when one allows for any A' the perfomance 
of the algorithm is extremely poor: with the simple Gaussian 
measure for lattices ^ exp(— 7TrAj4*) we were able 

to recover the best packers in d = 2, 3, although already in 3 
dimensions we had to go to very low termperatures. The per- 



formance of the algorithm quickly deteriorates with dimen- 
sions, and by d = 10 it is completely useless. The above men- 
tioned variant of the algorithm where one samples only among 
the non-isometric lattices has similar performance when ap- 
proximantions are used. Finally it's worth mentioning that the 
lattices generated by such Markov chains are never perfect and 
are typically far from being such. 

These negative results provide an extra motivation for 
studying perfect lattices and the Voronoi construction where 
much better performance is achieved. 
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